#!/usr/bin/perl -w
use strict;
use Getopt::Long;
use lib "/net/cpp-group/Leo/bin";
use parse_bl2seq;
use parse_fasta;
#88888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Usage

#88888888888888888888888888888888888888888888888888888888888888888888888888888888
my $usage = <<'USAGE';

USAGE:

    fa_multiple_alignment_positions_with_gaps.pl
                            [--header]
                            [--help]
                            sequences.fa


        This script takes a multiple alignment
        prints the corresponding pairwise alignment results.
USAGE

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# mandatory


# optional parameters
my $help = undef;
my $header = undef;
GetOptions(
			'help' => \$help,
			'header' => \$header,

);

die $usage if ($help);

#_________________________________________________________________________________________

#   get_seq_callback
#                       $annotation     annotation line
#                       \@sequence      sequence lines
#                       $line count     line count in file
#                       $acc_pos        accession line file position
#                       $sequence_count index of current sequence
#                       $data           see below

#_________________________________________________________________________________________
my %sequences;
my $align_len;
sub get_seq_callback($$$$$$)
{
    my ($acc,  $seq, $save_sequences) = @_[0, 1, 5];

	my $sequence = join "", @$seq;
	$align_len = length($sequence);
	my @alignment = get_alignment_marks($sequence);
	$sequence =~ tr/\-//d;
	$sequences{$acc} = [gap_positions_array_to_str(@alignment), length($sequence)];
}

my $cnt_seq = 0;
if (!@ARGV)
{
	push(@ARGV, '-');
}
while (my $ARGV = shift @ARGV)
{
	unless (open (ARGV, $ARGV))
	{
		warn "Can't open $ARGV:$!\n";
		next;
	}
	$cnt_seq += parse_fasta::parse_sequences(  *ARGV,
												\&get_seq_callback,
												undef);
	close ARGV;
	
}

print STDERR "$cnt_seq\tsequences parsed.\n";
my @acc = keys %sequences;
if ($header)
{
	print
"seq1\tseq2\tstart\tlen\tpositions\tcoverage\te_value\tlen1\tlen2\tresidues\tbitscore\tgaps1\tgaps2\n".
"----\t----\t-----\t---\t---------\t--------\t-------\t----\t----\t--------\t--------\t-----\t-----\n";
}
for my $i (1.. $#acc)
{
	for my $j (0.. ($i-1))
	{
		my ($acc1, $acc2) = @acc[$i, $j];
		print join ("\t",
						$acc[$i],
						$acc[$j],
						0,
						$align_len,
						$align_len,
						1.0,
						0,
						$sequences{$acc1}[1],
						$sequences{$acc2}[1],
						$align_len,
						0,
						$sequences{$acc1}[0],
						$sequences{$acc2}[0]), "\n";
	}
}

